********************************************************************************
* Competition model solved using bi-level optimizaiton and branch & bound.
* Lei Zhang
* 10/28/2013
********************************************************************************

SETS
    i   Air_Products_facilities         /F1*F2/
    ii  Competitor_facilities           /FC1/
    j   Markets                         /M1*M2/
    k   Products                        /LOX,LNI/
    t   Time_periods                    /T1*T4/;

*-------------------------------------------------------------------------------

SCALAR
    CMAX    Upper_bound_of_capacity_expansion   /100.0/
    CMIN    Lower_bound_of_capacity_expansion   /5.0/
    MUU     Upper_bound_of_mu                   /1000.0/
    SU      Upper_bound_of_S                    /500.0/
    SCU     Upper_bound_of_SC                   /500.0/
    CCU     Upper_bound_of_CC                   /600.0/;

TABLE P(i,j,k,t)    Air_products_price 
            T1      T2      T3      T4
F1.M1.LOX   2.0     2.0     2.0     2.0
F1.M2.LOX   2.0     2.0     2.0     2.0
F1.M1.LNI   2.0     2.0     2.0     2.0
F1.M2.LNI   2.0     2.0     2.0     2.0
F2.M1.LOX   2.0     2.0     2.0     2.0
F2.M2.LOX   2.0     2.0     2.0     2.0
F2.M1.LNI   2.0     2.0     2.0     2.0
F2.M2.LNI   2.0     2.0     2.0     2.0;

TABLE PC(ii,j,k,t)  Competitor_price
            T1      T2      T3      T4
FC1.M1.LOX   2.1     2.1     2.1     2.1
FC1.M2.LOX   2.1     2.1     2.1     2.1
FC1.M1.LNI   2.1     2.1     2.1     2.1
FC1.M2.LNI   2.1     2.1     2.1     2.1;

TABLE FALPHA(i,k,t) Fixed_cost_to_expand_capacity
            T1      T2      T3      T4
F1.LOX      10.0    9.0     8.0     8.0
F1.LNI      10.0    9.0     8.0     8.0
F2.LOX      10.0    9.0     8.0     8.0
F2.LNI      10.0    9.0     8.0     8.0;

TABLE FCALPHA(ii,k,t)   Fixed_cost_for_competitor_to_expand_capacity
            T1      T2      T3      T4
FC1.LOX     10.1    9.1     8.1     8.1
FC1.LNI     10.1    9.1     8.1     8.1;

TABLE FBETA(i,k,t)  Unit_investment_cost_to_expand_capacity
            T1      T2      T3      T4
F1.LOX      1.1     0.9     0.8     0.8
F1.LNI      1.1     0.9     0.8     0.8
F2.LOX      1.1     0.9     0.8     0.8
F2.LNI      1.1     0.9     0.8     0.8;

TABLE FCBETA(ii,k,t)    Unit_investment_cost_for_competitor_to_expand_capacity
            T1      T2      T3      T4
FC1.LOX     1.2     1.0     0.9     0.8
FC1.LNI     1.2     1.0     0.9     0.8;

TABLE FGAMMA(i,k,t) Unit_production_cost
            T1      T2      T3      T4
F1.LOX      0.6     0.6     0.6     0.6
F1.LNI      0.6     0.6     0.6     0.6
F2.LOX      0.6     0.6     0.6     0.6
F2.LNI      0.6     0.6     0.6     0.6;

TABLE FCGAMMA(ii,k,t)   Unit_production_cost_for_competitor
            T1      T2      T3      T4
FC1.LOX     0.5     0.5     0.5     0.5
FC1.LNI     0.5     0.5     0.5     0.5;

TABLE TR(i,j,k) Unit_transportation_cost
            LOX     LNI
F1.M1       0.1     0.1
F1.M2       0.2     0.2
F2.M1       0.2     0.1
F2.M2       0.1     0.1;

TABLE TRC(ii,j,k)   Unit_transportation_cost_for_competitor
            LOX     LNI
FC1.M1      0.2     0.2
FC1.M2      0.1     0.1;

TABLE D(j,k,t)  Demand_of_market
            T1      T2      T3      T4
M1.LOX      35.0    38.0    40.0    42.0
M2.LOX      20.0    25.0    30.0    35.0
M1.LNI      35.0    38.0    40.0    42.0
M2.LNI      20.0    25.0    30.0    35.0;

PARAMETER
    xc(ii,k,t)
    temp
    temp1
    iter
    iter1
    iter2
    iter3
    n;

*-------------------------------------------------------------------------------

VARIABLES
    npv             Net_present_value
    npvc            Net_present_value_of_competitor
    lambda1(j,k,t)  KKT;

BINARY VARIABLES
    x(i,k,t)        Selection_of_capacity_expansion;

POSITIVE VARIABLES
    y(i,j,k,t)      Amount_of_product_that_Air_Products_sells_to_market
    yc(ii,j,k,t)    Amount_of_product_that_competitor_sells_to_market
    c(i,k,t)        Capacity_of_Air_products_facility
    dc(i,k,t)       Capacity_expansion_of_Air_products
    cc(ii,k,t)      Capacity_of_competitor_facility
    dcc(ii,k,t)     Capacity_expansion_of_competitor
    mu1(i,k,t)      KKT
    mu2(ii,k,t)     KKT
    s(i,k,t)        Total_amount_of_product_AP_sells_to_market
    sc(ii,k,t)      Total_amount_of_product_competitor_sells_to_market
    s1(i,k,t)       Convex_hull_of_s
    s2(i,k,t)       Convex_hull_of_s
    z1(i,k,t)       Convex_hull
    z2(ii,k,t)      Convex_hull
    sc1(ii,k,t)     Convex_hull_of_sc
    cc1(ii,k,t)     Convex_hull_of_cc
    cc2(ii,k,t)     Convex_hull_of_cc;

*-------------------------------------------------------------------------------

EQUATIONS
    objout          Objective_function_maximize_Air_Products_NPV
    objin           Objective_function_maximize_competitor_NPV
    epd(i,k,t)      Investment_decision_in_Air_Products_capacity_expansion
    upb(i,k,t)      Upper_bound_of_capacity_expansion
    epdc(ii,k,t)    Investment_decision_in_competitor_capacity_expansion
    stn1(i,j,k,t)   KKT_stationary
    stn2(ii,j,k,t)  KKT_stationary
    tmd(j,k,t)      Demand_satisfaction_for_all_markets
    spl(i,k,t)      Supply_of_Air_Products
    splc(ii,k,t)    Supply_of_competitor
    sch(i,k,t)      Supply_convex_hull
    dj1(i,k,t)      Disjunction1
    schc(ii,k,t)    Supply_competitor_convex_hull
    cchc(ii,k,t)    Capacity_competitor_convex_hull
    upbc(ii,k,t)    Upper_bound_of_capacity_expansion
    cs(i,k,t)       Supply_and_capacity
    csc(ii,k,t)     Supply_and_competitor_capacity
    mub(i,k,t)      Upper_bound_mu
    sub(i,k,t)      Upper_bound_s
    sub2(i,k,t)     Upper_bound_s2
    mub2(ii,k,t)    Upper_bound_mu2
    scub(ii,k,t)    Upper_bound_sc1
    scub2(ii,k,t)   Upper_bound_sc2
    ccub(ii,k,t)    Upper_bound_cc1
    ccub2(ii,k,t)   Upper_bound_cc2
    zub1(i,k,t)     Upper_bound_z1
    zub2(ii,k,t)    Upper_bound_z2
    ;


objout..            npv =e= SUM((i,j,k,t), P(i,j,k,t)*y(i,j,k,t))
                        - SUM((i,k,t), FALPHA(i,k,t)*x(i,k,t))
                        - SUM((i,k,t), FBETA(i,k,t)*dc(i,k,t))
                        - SUM((i,k,t), FGAMMA(i,k,t)*c(i,k,t))
                        - SUM((i,j,k,t), TR(i,j,k)*y(i,j,k,t));
objin..             npvc =e= SUM((ii,j,k,t), PC(ii,j,k,t)*yc(ii,j,k,t))
                        - SUM((ii,k,t), FCALPHA(ii,k,t)*xc(ii,k,t))
                        - SUM((ii,k,t), FCBETA(ii,k,t)*dcc(ii,k,t))
                        - SUM((ii,k,t), FCGAMMA(ii,k,t)*cc(ii,k,t))
                        - SUM((ii,j,k,t), TRC(ii,j,k)*yc(ii,j,k,t));

epd(i,k,t)..        c(i,k,t) =e= c(i,k,t-1) + dc(i,k,t);
upb(i,k,t)..        dc(i,k,t) =l= x(i,k,t) * CMAX;
epdc(ii,k,t)..      cc(ii,k,t) =e= cc(ii,k,t-1) + dcc(ii,k,t);
stn1(i,j,k,t)..     P(i,j,k,t) + mu1(i,k,t) + lambda1(j,k,t) =e= 0.0;
stn2(ii,j,k,t)..    PC(ii,j,k,t) + mu2(ii,k,t) + lambda1(j,k,t) =e= 0.0;
tmd(j,k,t)..        SUM(i, y(i,j,k,t)) + SUM(ii, yc(ii,j,k,t)) =e= D(j,k,t);
spl(i,k,t)..        s(i,k,t) =e= SUM(j, y(i,j,k,t));
splc(ii,k,t)..      sc(ii,k,t) =e= SUM(j, yc(ii,j,k,t));
sch(i,k,t)..        s(i,k,t) =e= s1(i,k,t) + s2(i,k,t);
dj1(i,k,t)..        s2(i,k,t) =e= c(i,k,t) * (1.0-z1(i,k,t));
schc(ii,k,t)..      sc(ii,k,t) =e= sc1(ii,k,t) + cc2(ii,k,t);
cchc(ii,k,t)..      cc(ii,k,t) =e= cc1(ii,k,t) + cc2(ii,k,t);
upbc(ii,k,t)..      dcc(ii,k,t) =l= CMAX * xc(ii,k,t);
cs(i,k,t)..         s(i,k,t) =l= c(i,k,t);
csc(ii,k,t)..       sc(ii,k,t) =l= cc(ii,k,t);
mub(i,k,t)..        mu1(i,k,t) =l= MUU * (1.0-z1(i,k,t));
sub(i,k,t)..        s1(i,k,t) =l= SU * z1(i,k,t);
sub2(i,k,t)..       s2(i,k,t) =l= SU * (1.0-z1(i,k,t));
mub2(ii,k,t)..      mu2(ii,k,t) =l= MUU * (1.0-z2(ii,k,t));
scub(ii,k,t)..      sc1(ii,k,t) =l= SCU * z2(ii,k,t);
scub2(ii,k,t)..     cc2(ii,k,t) =l= SCU * (1.0-z2(ii,k,t));
ccub(ii,k,t)..      cc1(ii,k,t) =l= CCU * z2(ii,k,t);
ccub2(ii,k,t)..     cc2(ii,k,t) =l= CCU * (1.0-z2(ii,k,t));
zub1(i,k,t)..       z1(i,k,t) =l= 1.0;
zub2(ii,k,t)..      z2(ii,k,t) =l= 1.0;

MODEL COMPETITION /all/;

$echo bilevel c dc x max npvc objin epdc stn1 stn2 tmd spl splc sch dj1 schc cchc upbc cs csc mub sub sub2 mub2 scub scub2 ccub ccub2 zub1 zub2 > "%emp.info%";

xc(ii,k,t) = 0;

temp1 = -100000;
n = 0;

file results /results.dat/;
put results;
put '        Optimization resuls of tri-level competition problem.'//;
put '                                          date: ',system.date/;
put '                                          time: ',system.time/;
put /;

put '        No. of our facilities: ', card(i)/;
put '        No. of competitor facilities: ', card(ii)/;
put '        No. of products: ', card(k)/;
put '        No. of time periods: ', card(t);

put //;
put '          No.'
loop((ii,k,t), put '          xc');
put '         npv';
put '        npvc';
loop((i,k,t), put '           c');
loop((ii,k,t), put '          cc');
put /;

for(iter = 1 to power(2,card(ii)*card(k)*card(t)),
    temp = iter;
    for(iter1 = 1 to card(ii),
        for(iter2 = 1 to card(k),
            for(iter3 = 1 to card(t),
                xc(ii,k,t)$((ord(ii)=iter1) and (ord(k)=iter2) and (ord(t)=iter3))
                    = mod(temp,2);
                temp = floor(temp/2);
                );
            );
        );
    SOLVE COMPETITION USING emp MAXIMIZING npv;
    if(temp1 le npvc.l,
        temp1 = npvc.l;
        n = iter;
        );
    put iter;
    loop((ii,k,t), put xc(ii,k,t));
    put npv.l;
    put npvc.l;
    loop((i,k,t), put c.l(i,k,t));
    loop((ii,k,t), put cc.l(ii,k,t));
    put /;
    );

put /;
put '*** No. ', n, ' is the optimal solution.'
put //;

